﻿module ExpCount

    let benchmark f loops = 
        let timer = new System.Diagnostics.Stopwatch()
        timer.Start()
        let returnValue =
            for i = 0 to loops do f()
        printfn "Elapsed Time: %i" timer.ElapsedMilliseconds
    let isqrt =
        let rec f i c n =
            match n<i with
            | false -> f (i+2) (c+1) (n-i) 
            | true -> c
        f 1 0
    let IssueError msg = printfn "%s" msg
    let unknown = '?'

    // Caculates the hamming distance of the known values between two char arrays, <a> and <b>
    let HammingDistance (a : char array) (b : char array) = 
        match (a, b) with
        | a, b when a.Length <> b.Length -> IssueError "The two char arrays are of different lengths"; -1
        | t -> t ||> Seq.map2 (fun i j -> if i = unknown || j = unknown || i = j then 0 else 1) |> Seq.sum

    (* let ShiftingHammingDistance a b (k : int) = 
        if a.Length <> b.Length then
            -1
        else
            let length = a.Length
            let mutable count = 0

            for i in length - k .. length do
                //printf "%s - %s" (a.Substring(0, i)) (b.Substring(length - i))
                //printf ":%d\n" (HammingDistance (a.Substring(0, i)) (b.Substring(length - i)))


                // CHECK IF EVERYTHING REMAINS IN BOUNDS
                if (HammingDistance (a[0 .. i]) (b.[(length - i) .. length])) >= i - length + k then
                    count <- count + 1
            count
        


    let factorial a =
        //| a when a <= 1 -> 1
        //| _ -> a * factorial (a - 1)
        *)

    // A choose b
    let combinations a b =
        match (a, b) with
        | (a, b) when b > a || b < 0 -> 0
        | _ -> [| 0 .. b - 1 |] |> Array.fold (fun acc i -> acc * (a - i) / (i + 1)) 1
            
        

    // Exp1 [|'a';'a';'?';'?'|] [|'a';'?';'a';'?'|] 2 4;;

    // Calculates the p probability that <a> and <b> are meet constraint 1 of <k>
    //  - Constraint 1 is the minimum hamming distance, p is calculated as 1/gluenumber
    let Probability1 (a : char array) (b : char array) (k : int) (glueNumber : int) =
        match a, b, k, glueNumber with
        | a, b, _, _ when a.Length <> b.Length -> IssueError "The two char arrays are of different lengths"; -1.0
        | _, _, k, _ when k < 0 -> IssueError "k, the constraint, cannot be negative"; -1.0
        | _, _, _, g when g <= 0 -> IssueError "The number of different glues must be greater than 0"; -1.0
        | _ ->
            let distanceLeft = k - (HammingDistance a b) // counts the hamming distance disregarding unknowns
            // Progresses along <a> and <b> by pairs and counting the number of those pairs that contain an unknown in either side
            let unknownPairs = // e.g. "a?a?" "aa??"would be 3 if 'q' is the value of unknown
                [| 0 .. a.Length - 1 |]
                |> Array.sumBy (fun i -> if (a.[i] = unknown || b.[i] = unknown) then 1 else 0)
            let p = (1.0 / (float glueNumber)) // The probability of it not occuring
            if distanceLeft > 0 then
                let c =
                    // We start at distanceLeft because it will save memory depending on how combinations is done
                    // This is so that we only created as much memory as needed
                    (([|combinations unknownPairs distanceLeft|], [| distanceLeft .. unknownPairs |])
                    ||> Array.fold (fun acc i ->
                        let last = Array.get acc (i - distanceLeft)
                        Array.append acc [| last * (unknownPairs - i) / (i + 1)|]))
                
                // Similar to the binary formula, summing over from distanceleft to unknownPairs:
                // (unknown pairs choose i) * p^i * pbar^(unknownPairs - i)
                //TODO: Check the accuracy of this lower bound and pbar
                (0.0, [| distanceLeft .. unknownPairs|]) ||> Array.fold (fun acc i ->
                    acc + (float c.[i - distanceLeft]) * (p  ** (float i)) * ((1.0 - p) ** (float <| unknownPairs - i)))
            else // Means that it already matches the 
                1.0

    // Calculates the probability that <a> and <b> are meet constraint 4 of <k> at the <i> index
    //  - Constraint 4 is the minimum shifting hamming distance, p is calculated as 1/gluenumber
    let Probability4 (a : char array) (b : char array) (k : int) (glueNumber : int) (i : int) =
        match (a, b, k, glueNumber, i) with
        | a, b, _, _, _ when a.Length <> b.Length -> IssueError "The two char arrays are of different lengths"; -1.0
        | _, _, k, _, _ when k < 0 -> IssueError "k, the constraint, cannot be negative"; -1.0
        | _, _, _, g, _ when g <= 0 -> IssueError "The number of different glues must be greater than 0"; -1.0
        | a, b, k, _, i when i < (a.Length - 1 - k) || i >= a.Length -> -1.0 //IssueError "incorrect i index"; -1.0
        | _ -> Probability1 a.[0 .. i] b. [a.Length - 1 - i .. a.Length - 1] k glueNumber
        
    (*
    let bin = System.Random()
    let test = Array2D.init 5 8 (fun i j ->  [|'a'; 'b'; 'c'; 'd';'?'|].[bin.Next(5)])
    ExpCount test 2 4 4
    *)
    (*
    let mutable testM = Array2D.init 2 5 (fun i j -> unknown)
    testM.[0,0] <- 'a'
    testM.[1,0] <- 'a'
    testM.[0,1] <- 'b'
    testM.[1,1] <- 'a'
    testM.[0,2] <- 'c'
    testM.[1,2] <- 'a'
    testM.[0,3] <- 'd'
    testM.[1,3] <- 'a'
    testM.[0,4] <- 'a'
    testM.[1,4] <- 'a'
    ExpCount testM 1 1 4
    *)

    /// Returns the probabilty that two strings a and b, with given unknowns will satisfy the two contraints
    /// C1 (Hammond Distance >= k1)  and C4 (Shifting Hammond Distance >= k4)
    let ExpCount (M : char [,]) k1 k4 glueNumber =
        let k = max k1 k4
        let el = M.GetLength 0
        let n = M.GetLength 1
        (*let findIndex (i : int) =
            let x = (-1 + (int <| sqrt (1.0 + 8.0 * (float i)))) / 2
            //let x = (-1 + (isqrt 1 + 8 * i)) / 2
            (el - 2 - x, el - 1 - i + (x * (x + 1) / 2))
        // findIndex is designed to replace the following:
        [| 0 .. (el - 1) * el / 2 - 1 |]
            |> Array.map (fun y ->
                let (i, j) = findIndex y*)
        [| 0 .. n - 2 |] |> Array.map (fun i ->
            [| i + 1 .. n - 1 |] |> Array.map (fun j ->
                let Mi = Array.init el (fun x -> M.[x, i])
                let Mj = Array.init el (fun x -> M.[x, j])
                let e2 =
                    [| n - k4 .. n - 2 |]
                    |> Array.map(fun x -> 2.0 - (Probability4 Mi Mj k4 glueNumber x) + (Probability4 Mj Mi k4 glueNumber x))
                    |> Array.sum
                let e1 = (Probability1 Mi Mj k glueNumber)
                //printfn "%f: %A<->%A" e1 Mi Mj
                e1 + e2
                ) |> Array.sum)
            |> Array.sum

    (*
    let test a = 
        DetWords 5 10 3 2 [|'a';'b';'c';'d'|]
        printf ""
    benchmark test 1
    *)

    // No gareentee that any of the glues will work
    // Specify <glues> in order of preferance from most preferred to least preferred
    // - <el> is the length of the words
    // - <n> is the number of words to be produced
    // - <k1> is the C1 constraint, the minimum hamming distance between all the produced dna words
    // - <k4> is the C4 constraint, the minimum shifting hamming distance between all the produced dna words
    // - <glues> is the array of possible entites of which a dna word can consist
    let DetWords n el k1 k4 (glues : char array) =
        let mutable M = Array2D.init el n (fun i j -> unknown)
        
        for y = 0 to n - 1 do
            for x = 0 to el - 1 do
                
                // Finds the glue with the largest expectation
             //   if (ExpCount M k1 k4 glues.Length) = 0.0 then
             //       IssueError "No possible word set will fit the constraints"
                printfn "%d, %d:" x y
                let mutable strongestGlue = unknown
                let mutable expectation = 0.0
                for i = 0 to glues.Length - 1 do
                    M.[x, y] <- glues.[i]
                    //let e = ExpCount M.[0 .. el - 1, 0 .. y + if y < n - 1 then 1 else 0] k1 k4 glues.Length
                    let e = ExpCount M k1 k4 glues.Length
                    printfn "%c: %f" glues.[i] e
                    if e > expectation then
                        strongestGlue <- glues.[i]
                        expectation <- e
                M.[x, y] <- strongestGlue
        M




    // Section 3.3, pg 14, lemma 10
    // Calculates the float value for an expectation as illustrated in
    let ExpCount2 n el k1 k4 =
        if el < k4 then
            IssueError "el must be greater than or equal to k4"
            -1.0
        elif el <= 0 then
            IssueError "el must be greater than 0"
            -1.0
        elif n < 0 then
            IssueError "n cannot be negative"
            -1.0
        elif (k1 < 0) || (k4 < 0) then
            IssueError "k1 and k4, the constraints, cannot be negative"
            -1.0
        else
            let mutable sum1 = 0.0
            let mutable c = 1 // will hold the combinations and j starts at 0, anything choose 0 is 1 and el is not zero
            for j = 0 to (max k1 k4) - 1 do
                sum1 <- sum1 + (float <| c) * (2.0 ** (float -el))
                c <- c * (el - j) / (j + 1)
        
            let mutable sum2 = 0.0
            for i = el - k4 + 1 to el - 1 do
                c <- 1 // j starts at 0, anything choose 0 is 1 and i is never 0
                for j = 0 to k4 - el + i - 1 do
                    sum2 <- sum2 + (float <| c) * (2.0 ** (float -i))
                    c <- c * (i - j) / (j + 1)
            (sum1 + 2.0 * sum2) * (float <| combinations n 2)

    // general function for binomial coefficients
    // TODO: do the division separately
    let choose n k = List.fold (fun s i -> s * (n-i+1)/i ) 1 [1..k]
    //a recursive binary search for optimal l
    let rec findBinaryL (l_star:int) (n:int) (k1:int) (k4:int) =
        let k = max k1 k4
        if l_star < k then
            -1 // return bad value
        elif n <= 0 then
            -1
        else
            let l_separator = k + ((l_star-k)/2) // l_star is the initial upper bound. we will feed this constantly where k is the minimum
            let right = float((choose n 2)*(1+2*(k4-1)) - 1)
            let left = ExpCount2 n l_separator k1 k4
            if left > right then // l_star can be minimized further, so feed l_separator as the upper bound
                findBinaryL l_separator n k1 k4
            elif (l_star-l_separator) > 1 then // is this the right pre-termination step?
                findBinaryL l_star n k1 l_separator
            else
                l_separator


    // Section 3.3, pg 14, lemma 11
    // Runs a binary search looking the for optimal L for a binary glue set
    let findOptimalL n l_star k1 k4 = 
        if l_star < k4 then
            IssueError "l_star must be greater than or equal to k4"
            -1
        elif l_star <= 0 then
            IssueError "l_star must be greater than 0"
            -1
        elif n < 0 then
            IssueError "n cannot be negative"
            -1
        elif (k1 < 0) || (k4 < 0) then
            IssueError "k1 and k4, the constraints, cannot be negative"
            -1
        else
            let mutable lowerBound = k4
            let mutable upperBound = l_star
            let mutable trial = (lowerBound + upperBound + 1) / 2    // plus one for integer ciel
            let mutable lastTrial = 0
            //let testValue = float <| (1 + 2 * (k4 - 1)) * (combinations n 2) - 1
            let testValue = 1.0
        
            let epsilon = 0.05

            //printfn "OptimalL %f:%f" testValue <| ExpCount2 n trial k1 k4

            // Run the loop that performs the binary search until the last guest is the same
            // The commented section is a test to see if the difference between the two floats is insignificant
            while (trial <> lastTrial) (*&& (abs <| (ExpCount2 n trial k1 k4 ) - testValue > epsilon) *) do
                lastTrial <- trial

                //printfn "OptimalL %f:%f" testValue <| ExpCount2 n trial k1 k4
            
                // If the guess produces a larger value than the test, guess smaller, else guess larger
                if (ExpCount2 n trial k1 k4) < testValue then
                    upperBound <- trial
                    trial <- (lowerBound + upperBound + 1) / 2   // plus one for integer ciel
                    //printfn "a:%d:%d" lowerBound upperBound
                else
                    lowerBound <- trial
                    trial <- (lowerBound + upperBound + 1) / 2   // plus one for integer ciel
                    //printfn "b:%d:%d" lowerBound upperBound
            trial

    // Section 3.3, pg 14, lemma 11
    // Runs a binary search looking the for optimal L for a binary glue set
    let findExhuastiveOptimalL n l_star k1 k4 = 
        if l_star < k4 then
            IssueError "l_star must be greater than or equal to k4"
            -1
        elif l_star <= 0 then
            IssueError "l_star must be greater than 0"
            -1
        elif n < 0 then
            IssueError "n cannot be negative"
            -1
        elif (k1 < 0) || (k4 < 0) then
            IssueError "k1 and k4, the constraints, cannot be negative"
            -1
        else
            let mutable trial = k4
            let mutable expectation = ExpCount2 n trial k1 k4
            let testValue = 1.0 //float <| -1 + (2 * k4 - 1) * (combinations n 2)
        
            printfn "OptimalL %d:%f" trial <| ExpCount2 n trial k1 k4
            while (ExpCount2 n trial k1 k4) >= testValue do
                trial <- trial + 1
                printfn "OptimalL %d:%f" trial <| ExpCount2 n trial k1 k4
            trial

(* Figure out E_1, E_4 *)
(* E_1 is the count of C1 bool true; E_4 is the count of C4 bool true; make histograms

E4 = related to distance

distance 4  C(3 2) * (1/2)^3
distance 5  C(3 3) * (1/2)^3
*)